Setup

## Directory setup
here_root <- "benchmarks/lupus-n_patients"
here::i_am(file.path(here_root, "analysis/lupus-n_patients-mock-results.Rmd"))
#> here() starts at /kyukon/data/gent/vo/000/gvo00063/Jeroen/DD_project_v2

res_dir <- here::here(here_root, "results")
fig_dir <- here::here(here_root, "figures")

Using n_patients: 10 20 30
Using methods: edgeR_NB, edgeR_QP, edgeR_NB_optim, edgeR_QP_optim, bGLM, qbGLM, qbGLM_offset, qbGLM_offset_squeeze
Using celltype: B_mem

Load results

# TODO: also add all-patients results for comparison?



res_files <- map(n_patients, ~ get_mock_res_files(
    dataset = "lupus-n_patients",
    methods = methods, celltype = celltype, datatype = "pb",
    n_patients = .x
)) %>%
    set_names(paste0("n_patients_", n_patients))

#res_list <- map_depth(res_files, 2, readRDS)
res_list <- lapply(res_files, function(element){
    lapply(element, readRDS)
})
str(res_list, max.level = 3)
#> List of 3
#>  $ n_patients_10:List of 8
#>   ..$ edgeR_NB            :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ edgeR_QP            :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ edgeR_NB_optim      :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ edgeR_QP_optim      :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ bGLM                :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ qbGLM               :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ qbGLM_offset        :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ qbGLM_offset_squeeze:List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>  $ n_patients_20:List of 8
#>   ..$ edgeR_NB            :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ edgeR_QP            :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ edgeR_NB_optim      :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ edgeR_QP_optim      :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ bGLM                :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ qbGLM               :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ qbGLM_offset        :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ qbGLM_offset_squeeze:List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>  $ n_patients_30:List of 8
#>   ..$ edgeR_NB            :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ edgeR_QP            :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ edgeR_NB_optim      :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ edgeR_QP_optim      :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ bGLM                :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ qbGLM               :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ qbGLM_offset        :List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2
#>   ..$ qbGLM_offset_squeeze:List of 5
#>   .. ..$ replicate_1:List of 2
#>   .. ..$ replicate_2:List of 2
#>   .. ..$ replicate_3:List of 2
#>   .. ..$ replicate_4:List of 2
#>   .. ..$ replicate_5:List of 2

Load SCE objects

data_files <- map(n_patients, ~ get_SCE_files(
    dataset = "lupus-n_patients", which = "mock_replicates",
    celltype = celltype, n_patients = .x
)) %>%
    set_names(paste0("n_patients_", n_patients))
sce_objects <- map(data_files, readRDS)
str(sce_objects, max.level = 3)
#> List of 3
#>  $ n_patients_10:List of 5
#>   ..$ replicate_1:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_2:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_3:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_4:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_5:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>  $ n_patients_20:List of 5
#>   ..$ replicate_1:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_2:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_3:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_4:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_5:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>  $ n_patients_30:List of 5
#>   ..$ replicate_1:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_2:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_3:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_4:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
#>   ..$ replicate_5:Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots

Data overview

The results were generated from 5 mock replicates
Each replicate was generated by randomly splitting the subjects in two mock groups
No sub-sampling of cells per patient was performed for this data
The number of patients was randomly sub-sampled stratified by mock group

Note that for n_patients 10 and 20, the selected patients vary randomly between the replicates, which is why the nrows and ncols are not identical. For the n_patients 30, however, no random sub-sampling was done because I could just use all patients from 2 out of the original 3 batches from the Control samples. That’s why the nrows and ncols are identical for those subsets.

The mock group assignment of each subject also varies across replicates.

map_dfr(sce_objects,
    ~ map_dfr(.x, function(x) c(nrows = nrow(x), ncols = ncol(x)),
        .id = "replicate"
    ),
    .id = "n_patients"
)

Subjects are divided across mock groups as follows:

`n_patients = 10`

replicate_1

	A	B
IGTB1762_IGTB1762	1	0
IGTB1793_IGTB1793	0	1
IGTB1815_IGTB1815	0	1
IGTB1828_IGTB1828	1	0
IGTB1871_IGTB1871	0	1
IGTB1895_IGTB1895	1	0
IGTB1982_IGTB1982	0	1
IGTB1996_IGTB1996	1	0
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	0	1

replicate_2

	A	B
IGTB1762_IGTB1762	0	1
IGTB1789_IGTB1789	0	1
IGTB1828_IGTB1828	1	0
IGTB1840_IGTB1840	0	1
IGTB1871_IGTB1871	1	0
IGTB1895_IGTB1895	0	1
IGTB1901_IGTB1901	1	0
IGTB1952_IGTB1952	1	0
IGTB1996_IGTB1996	0	1
IGTB2007_IGTB2007	1	0

replicate_3

	A	B
IGTB1762_IGTB1762	0	1
IGTB1789_IGTB1789	1	0
IGTB1793_IGTB1793	1	0
IGTB1815_IGTB1815	0	1
IGTB1895_IGTB1895	1	0
IGTB1952_IGTB1952	0	1
IGTB1966_IGTB1966	1	0
IGTB1996_IGTB1996	0	1
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	0	1

replicate_4

	A	B
IGTB1762_IGTB1762	0	1
IGTB1789_IGTB1789	1	0
IGTB1815_IGTB1815	1	0
IGTB1828_IGTB1828	0	1
IGTB1895_IGTB1895	0	1
IGTB1901_IGTB1901	1	0
IGTB1952_IGTB1952	0	1
IGTB1982_IGTB1982	0	1
IGTB1996_IGTB1996	1	0
IGTB2065_IGTB2065	1	0

replicate_5

	A	B
IGTB1762_IGTB1762	1	0
IGTB1789_IGTB1789	1	0
IGTB1793_IGTB1793	0	1
IGTB1815_IGTB1815	0	1
IGTB1828_IGTB1828	0	1
IGTB1840_IGTB1840	0	1
IGTB1901_IGTB1901	1	0
IGTB1996_IGTB1996	0	1
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	1	0

`n_patients = 20`

replicate_1

	A	B
IGTB141_IGTB141	1	0
IGTB143_IGTB143	1	0
IGTB195_IGTB195	0	1
IGTB469_IGTB469	1	0
IGTB498_IGTB498	1	0
IGTB508_IGTB508	1	0
IGTB514_IGTB514	1	0
IGTB826_IGTB826	0	1
IGTB1506_IGTB1506	0	1
IGTB1762_IGTB1762	1	0
IGTB1789_IGTB1789	0	1
IGTB1793_IGTB1793	0	1
IGTB1840_IGTB1840	0	1
IGTB1871_IGTB1871	0	1
IGTB1895_IGTB1895	1	0
IGTB1952_IGTB1952	1	0
IGTB1966_IGTB1966	0	1
IGTB1982_IGTB1982	0	1
IGTB1996_IGTB1996	1	0
IGTB2065_IGTB2065	0	1

replicate_2

	A	B
IGTB141_IGTB141	1	0
IGTB143_IGTB143	0	1
IGTB195_IGTB195	1	0
IGTB469_IGTB469	1	0
IGTB508_IGTB508	0	1
IGTB514_IGTB514	0	1
IGTB670_IGTB670	0	1
IGTB826_IGTB826	1	0
IGTB1506_IGTB1506	0	1
IGTB1575_IGTB1575	1	0
IGTB1762_IGTB1762	0	1
IGTB1789_IGTB1789	0	1
IGTB1793_IGTB1793	0	1
IGTB1840_IGTB1840	0	1
IGTB1871_IGTB1871	1	0
IGTB1901_IGTB1901	1	0
IGTB1966_IGTB1966	0	1
IGTB1982_IGTB1982	1	0
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	1	0

replicate_3

	A	B
IGTB143_IGTB143	0	1
IGTB195_IGTB195	0	1
IGTB498_IGTB498	0	1
IGTB508_IGTB508	0	1
IGTB514_IGTB514	0	1
IGTB645_IGTB645	0	1
IGTB670_IGTB670	1	0
IGTB826_IGTB826	1	0
IGTB1372_IGTB1372	1	0
IGTB1539_IGTB1539	1	0
IGTB1575_IGTB1575	1	0
IGTB1789_IGTB1789	1	0
IGTB1815_IGTB1815	0	1
IGTB1828_IGTB1828	1	0
IGTB1901_IGTB1901	1	0
IGTB1906_IGTB1906	1	0
IGTB1952_IGTB1952	0	1
IGTB1982_IGTB1982	0	1
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	0	1

replicate_4

	A	B
IGTB143_IGTB143	1	0
IGTB469_IGTB469	0	1
IGTB498_IGTB498	0	1
IGTB670_IGTB670	1	0
IGTB1372_IGTB1372	0	1
IGTB1506_IGTB1506	0	1
IGTB1539_IGTB1539	1	0
IGTB1575_IGTB1575	1	0
IGTB1762_IGTB1762	0	1
IGTB1789_IGTB1789	1	0
IGTB1840_IGTB1840	0	1
IGTB1871_IGTB1871	0	1
IGTB1895_IGTB1895	0	1
IGTB1901_IGTB1901	1	0
IGTB1906_IGTB1906	1	0
IGTB1952_IGTB1952	0	1
IGTB1966_IGTB1966	1	0
IGTB1982_IGTB1982	0	1
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	1	0

replicate_5

	A	B
IGTB141_IGTB141	0	1
IGTB143_IGTB143	0	1
IGTB498_IGTB498	1	0
IGTB508_IGTB508	0	1
IGTB514_IGTB514	1	0
IGTB645_IGTB645	0	1
IGTB1506_IGTB1506	0	1
IGTB1539_IGTB1539	1	0
IGTB1575_IGTB1575	1	0
IGTB1762_IGTB1762	1	0
IGTB1789_IGTB1789	1	0
IGTB1815_IGTB1815	0	1
IGTB1840_IGTB1840	0	1
IGTB1895_IGTB1895	1	0
IGTB1901_IGTB1901	1	0
IGTB1906_IGTB1906	0	1
IGTB1966_IGTB1966	1	0
IGTB1982_IGTB1982	0	1
IGTB1996_IGTB1996	0	1
IGTB2065_IGTB2065	1	0

`n_patients = 30`

replicate_1

	A	B
IGTB141_IGTB141	1	0
IGTB143_IGTB143	1	0
IGTB195_IGTB195	0	1
IGTB469_IGTB469	1	0
IGTB498_IGTB498	1	0
IGTB508_IGTB508	1	0
IGTB514_IGTB514	1	0
IGTB645_IGTB645	1	0
IGTB670_IGTB670	0	1
IGTB826_IGTB826	0	1
IGTB1372_IGTB1372	0	1
IGTB1506_IGTB1506	0	1
IGTB1539_IGTB1539	1	0
IGTB1575_IGTB1575	0	1
IGTB1762_IGTB1762	1	0
IGTB1789_IGTB1789	0	1
IGTB1793_IGTB1793	0	1
IGTB1815_IGTB1815	0	1
IGTB1828_IGTB1828	1	0
IGTB1840_IGTB1840	0	1
IGTB1871_IGTB1871	0	1
IGTB1895_IGTB1895	1	0
IGTB1901_IGTB1901	1	0
IGTB1906_IGTB1906	0	1
IGTB1952_IGTB1952	1	0
IGTB1966_IGTB1966	0	1
IGTB1982_IGTB1982	0	1
IGTB1996_IGTB1996	1	0
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	0	1

replicate_2

	A	B
IGTB141_IGTB141	1	0
IGTB143_IGTB143	0	1
IGTB195_IGTB195	1	0
IGTB469_IGTB469	1	0
IGTB498_IGTB498	0	1
IGTB508_IGTB508	0	1
IGTB514_IGTB514	0	1
IGTB645_IGTB645	0	1
IGTB670_IGTB670	0	1
IGTB826_IGTB826	1	0
IGTB1372_IGTB1372	0	1
IGTB1506_IGTB1506	0	1
IGTB1539_IGTB1539	1	0
IGTB1575_IGTB1575	1	0
IGTB1762_IGTB1762	0	1
IGTB1789_IGTB1789	0	1
IGTB1793_IGTB1793	0	1
IGTB1815_IGTB1815	1	0
IGTB1828_IGTB1828	1	0
IGTB1840_IGTB1840	0	1
IGTB1871_IGTB1871	1	0
IGTB1895_IGTB1895	0	1
IGTB1901_IGTB1901	1	0
IGTB1906_IGTB1906	1	0
IGTB1952_IGTB1952	1	0
IGTB1966_IGTB1966	0	1
IGTB1982_IGTB1982	1	0
IGTB1996_IGTB1996	0	1
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	1	0

replicate_3

	A	B
IGTB141_IGTB141	0	1
IGTB143_IGTB143	0	1
IGTB195_IGTB195	0	1
IGTB469_IGTB469	1	0
IGTB498_IGTB498	0	1
IGTB508_IGTB508	0	1
IGTB514_IGTB514	0	1
IGTB645_IGTB645	0	1
IGTB670_IGTB670	1	0
IGTB826_IGTB826	1	0
IGTB1372_IGTB1372	1	0
IGTB1506_IGTB1506	1	0
IGTB1539_IGTB1539	1	0
IGTB1575_IGTB1575	1	0
IGTB1762_IGTB1762	0	1
IGTB1789_IGTB1789	1	0
IGTB1793_IGTB1793	1	0
IGTB1815_IGTB1815	0	1
IGTB1828_IGTB1828	1	0
IGTB1840_IGTB1840	1	0
IGTB1871_IGTB1871	0	1
IGTB1895_IGTB1895	1	0
IGTB1901_IGTB1901	1	0
IGTB1906_IGTB1906	1	0
IGTB1952_IGTB1952	0	1
IGTB1966_IGTB1966	1	0
IGTB1982_IGTB1982	0	1
IGTB1996_IGTB1996	0	1
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	0	1

replicate_4

	A	B
IGTB141_IGTB141	0	1
IGTB143_IGTB143	1	0
IGTB195_IGTB195	1	0
IGTB469_IGTB469	0	1
IGTB498_IGTB498	0	1
IGTB508_IGTB508	0	1
IGTB514_IGTB514	1	0
IGTB645_IGTB645	0	1
IGTB670_IGTB670	1	0
IGTB826_IGTB826	0	1
IGTB1372_IGTB1372	0	1
IGTB1506_IGTB1506	0	1
IGTB1539_IGTB1539	1	0
IGTB1575_IGTB1575	1	0
IGTB1762_IGTB1762	0	1
IGTB1789_IGTB1789	1	0
IGTB1793_IGTB1793	0	1
IGTB1815_IGTB1815	1	0
IGTB1828_IGTB1828	0	1
IGTB1840_IGTB1840	0	1
IGTB1871_IGTB1871	0	1
IGTB1895_IGTB1895	0	1
IGTB1901_IGTB1901	1	0
IGTB1906_IGTB1906	1	0
IGTB1952_IGTB1952	0	1
IGTB1966_IGTB1966	1	0
IGTB1982_IGTB1982	0	1
IGTB1996_IGTB1996	1	0
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	1	0

replicate_5

	A	B
IGTB141_IGTB141	0	1
IGTB143_IGTB143	0	1
IGTB195_IGTB195	0	1
IGTB469_IGTB469	1	0
IGTB498_IGTB498	1	0
IGTB508_IGTB508	0	1
IGTB514_IGTB514	1	0
IGTB645_IGTB645	0	1
IGTB670_IGTB670	1	0
IGTB826_IGTB826	1	0
IGTB1372_IGTB1372	1	0
IGTB1506_IGTB1506	0	1
IGTB1539_IGTB1539	1	0
IGTB1575_IGTB1575	1	0
IGTB1762_IGTB1762	1	0
IGTB1789_IGTB1789	1	0
IGTB1793_IGTB1793	0	1
IGTB1815_IGTB1815	0	1
IGTB1828_IGTB1828	0	1
IGTB1840_IGTB1840	0	1
IGTB1871_IGTB1871	0	1
IGTB1895_IGTB1895	1	0
IGTB1901_IGTB1901	1	0
IGTB1906_IGTB1906	0	1
IGTB1952_IGTB1952	0	1
IGTB1966_IGTB1966	1	0
IGTB1982_IGTB1982	0	1
IGTB1996_IGTB1996	0	1
IGTB2007_IGTB2007	1	0
IGTB2065_IGTB2065	1	0

Extract results of interest

Runtimes

runtimes <- map_dfr(res_list,
    ~ map_dfr(.x, get_runtimes, depth = 1, .id = "method"),
    .id = "n_patients"
) %>%
    mutate(n_patients = as.numeric(sub("n_patients_", "", n_patients)))

P-values

res_tables <- map_depth(
    res_list, 2,
    ~ get_aggregated_rep_tables(.x, depth = 1)
) %>%
    map(combine_tables, .id = "method")

Visualize results

Run times

ggplot(runtimes, aes(n_patients, time, col = method)) +
    geom_jitter(width = 0.5, alpha = 0.6) +
    geom_smooth(se = FALSE, method = "lm", formula = y ~ x) +
    labs(x = "Number of patients", y = "Time (seconds)", color = NULL) +
    scale_x_continuous(breaks = unique(runtimes$n_patients), minor_breaks = NULL)

Number of significant genes for different FDR levels

n_sign_summaries <- map(res_tables, function(tmp) {
    tmp %>%
        group_by(replicate, method) %>%
        summarise(
            sum(FDR < 0.01),
            sum(FDR < 0.05),
            sum(FDR < 0.1),
            .groups = "drop"
        )
})

n_patients_10

method	sum(FDR < 0.01)	sum(FDR < 0.05)	sum(FDR < 0.1)
Replicate: 1
bGLM	10	20	44
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	1
qbGLM_offset	1	1	5
qbGLM_offset_squeeze	0	0	0
Replicate: 2
bGLM	13	26	65
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 3
bGLM	12	31	68
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	1
edgeR_QP	0	0	0
edgeR_QP_optim	0	1	1
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 4
bGLM	30	72	109
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 5
bGLM	52	104	158
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0

n_patients_20

method	sum(FDR < 0.01)	sum(FDR < 0.05)	sum(FDR < 0.1)
Replicate: 1
bGLM	5	10	22
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 2
bGLM	23	57	85
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 3
bGLM	10	17	22
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 4
bGLM	9	12	14
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 5
bGLM	5	6	11
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0

n_patients_30

method	sum(FDR < 0.01)	sum(FDR < 0.05)	sum(FDR < 0.1)
Replicate: 1
bGLM	5	13	18
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 2
bGLM	27	55	81
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	1
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 3
bGLM	8	30	46
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 4
bGLM	11	21	25
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0
Replicate: 5
bGLM	8	13	18
edgeR_NB	0	0	0
edgeR_NB_optim	0	0	0
edgeR_QP	0	0	0
edgeR_QP_optim	0	0	0
qbGLM	0	0	0
qbGLM_offset	0	0	0
qbGLM_offset_squeeze	0	0	0

P-value distributions

pval_figs <- map(res_tables, pval_hist)

n_patients_10

n_patients_20

n_patients_30

P-value distributions stratified by expression level

For each replicate, calculate average log-normalized counts, using scuttle::calculateAverage().

## Calculate average counts for each replicate
avg_counts <- map_depth(sce_objects, 2,
    ~ calculateAverage(counts(.x))
)

## Split res_tables per replicate
res_tables_split <- map(res_tables, ~ split(.x, .x$replicate))

## Add average counts to results tables
## Note that for the 'combined' celltypes, the average counts will be the same
## for each celltype contrast
res_tables_split <- map2(
    res_tables_split, avg_counts,
    ~ map2(.x, .y, function(res, avg_cnt) {
        res$avg_count <- avg_cnt[res$gene]
        res
    })
)

Make p-value histograms stratified by the average counts.

pval_hist_strat_plots <- map(res_tables_split,
    ~ imap(.x, function(y, replicate) {
        replicate <- str_remove(replicate, "replicate_")
        title <- paste("Replicate:", replicate)
        pval_hist_strat(y, n_groups = 6) +
            ggtitle(title)
    })
)

n_patients_10

n_patients_20

n_patients_30

Session info

#> [1] "2023-06-02 20:24:26 CEST"
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.2.0 (2022-04-22)
#>  os       Red Hat Enterprise Linux 8.6 (Ootpa)
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Europe/Paris
#>  date     2023-06-02
#>  pandoc   2.13 @ /apps/gent/RHEL8/zen2-ib/software/Pandoc/2.13/bin/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  ! package              * version   date (UTC) lib source
#>  P argparse             * 2.2.2     2023-02-15 [?] CRAN (R 4.2.0)
#>  P beachmat               2.14.2    2023-04-07 [?] Bioconductor
#>  P beeswarm               0.4.0     2021-06-01 [?] CRAN (R 4.2.0)
#>  P Biobase              * 2.58.0    2022-11-01 [?] Bioconductor
#>  P BiocGenerics         * 0.44.0    2022-11-01 [?] Bioconductor
#>  P BiocManager            1.30.20   2023-02-24 [?] CRAN (R 4.2.0)
#>  P BiocNeighbors          1.16.0    2022-11-01 [?] Bioconductor
#>  P BiocParallel           1.32.6    2023-03-17 [?] Bioconductor
#>  P BiocSingular           1.14.0    2022-11-01 [?] Bioconductor
#>  P bitops                 1.0-7     2021-04-24 [?] CRAN (R 4.2.0)
#>  P bluster                1.8.0     2022-11-01 [?] Bioconductor
#>  P bslib                  0.4.2     2022-12-16 [?] CRAN (R 4.2.0)
#>  P cachem                 1.0.8     2023-05-01 [?] CRAN (R 4.2.0)
#>  P cli                    3.6.1     2023-03-23 [?] CRAN (R 4.2.0)
#>  P cluster                2.1.4     2022-08-22 [?] CRAN (R 4.2.0)
#>  P codetools              0.2-19    2023-02-01 [?] CRAN (R 4.2.0)
#>  P colorspace             2.1-0     2023-01-23 [?] CRAN (R 4.2.0)
#>    DDCompanion          * 0.1.15    2023-05-22 [1] local (./package)
#>  P DelayedArray           0.24.0    2022-11-01 [?] Bioconductor
#>  P DelayedMatrixStats     1.20.0    2022-11-01 [?] Bioconductor
#>  P digest                 0.6.31    2022-12-11 [?] CRAN (R 4.2.0)
#>  P dplyr                * 1.1.2     2023-04-20 [?] CRAN (R 4.2.0)
#>  P dqrng                  0.3.0     2021-05-01 [?] CRAN (R 4.2.0)
#>  P DT                     0.28      2023-05-18 [?] CRAN (R 4.2.0)
#>  P edgeR                  3.40.2    2023-01-19 [?] Bioconductor
#>  P ellipsis               0.3.2     2021-04-29 [?] CRAN (R 4.2.0)
#>  P evaluate               0.21      2023-05-05 [?] CRAN (R 4.2.0)
#>  P fansi                  1.0.4     2023-01-22 [?] CRAN (R 4.2.0)
#>  P farver                 2.1.1     2022-07-06 [?] CRAN (R 4.2.0)
#>  P fastmap                1.1.1     2023-02-24 [?] CRAN (R 4.2.0)
#>  P findpython             1.0.8     2023-03-14 [?] CRAN (R 4.2.0)
#>  P forcats              * 1.0.0     2023-01-29 [?] CRAN (R 4.2.0)
#>  P generics               0.1.3     2022-07-05 [?] CRAN (R 4.2.0)
#>  P GenomeInfoDb         * 1.34.9    2023-02-02 [?] Bioconductor
#>  P GenomeInfoDbData       1.2.9     2023-05-20 [?] Bioconductor
#>  P GenomicRanges        * 1.50.2    2022-12-16 [?] Bioconductor
#>  P ggbeeswarm             0.7.2     2023-04-29 [?] CRAN (R 4.2.0)
#>  P ggplot2              * 3.4.2     2023-04-03 [?] CRAN (R 4.2.0)
#>  P ggrepel                0.9.3     2023-02-03 [?] CRAN (R 4.2.0)
#>  P glue                   1.6.2     2022-02-24 [?] CRAN (R 4.2.0)
#>  P gridExtra              2.3       2017-09-09 [?] CRAN (R 4.2.0)
#>  P gtable                 0.3.3     2023-03-21 [?] CRAN (R 4.2.0)
#>  P here                 * 1.0.1     2020-12-13 [?] CRAN (R 4.2.0)
#>  P highr                  0.10      2022-12-22 [?] CRAN (R 4.2.0)
#>  P hms                    1.1.3     2023-03-21 [?] CRAN (R 4.2.0)
#>  P htmltools              0.5.5     2023-03-23 [?] CRAN (R 4.2.0)
#>  P htmlwidgets            1.6.2     2023-03-17 [?] CRAN (R 4.2.0)
#>  P httpuv                 1.6.11    2023-05-11 [?] CRAN (R 4.2.0)
#>  P httr                   1.4.6     2023-05-08 [?] CRAN (R 4.2.0)
#>  P iCOBRA                 1.26.0    2022-11-01 [?] Bioconductor
#>  P igraph                 1.4.3     2023-05-22 [?] CRAN (R 4.2.0)
#>  P IRanges              * 2.32.0    2022-11-01 [?] Bioconductor
#>  P irlba                  2.3.5.1   2022-10-03 [?] CRAN (R 4.2.0)
#>  P jquerylib              0.1.4     2021-04-26 [?] CRAN (R 4.2.0)
#>  P jsonlite               1.8.4     2022-12-06 [?] CRAN (R 4.2.0)
#>  P kableExtra           * 1.3.4     2021-02-20 [?] CRAN (R 4.2.0)
#>  P knitr                  1.42      2023-01-25 [?] CRAN (R 4.2.0)
#>  P labeling               0.4.2     2020-10-20 [?] CRAN (R 4.2.0)
#>  P later                  1.3.1     2023-05-02 [?] CRAN (R 4.2.0)
#>  P lattice                0.21-8    2023-04-05 [?] CRAN (R 4.2.0)
#>  P lifecycle              1.0.3     2022-10-07 [?] CRAN (R 4.2.0)
#>  P limma                  3.54.2    2023-02-28 [?] Bioconductor
#>  P locfit                 1.5-9.7   2023-01-02 [?] CRAN (R 4.2.0)
#>  P lubridate            * 1.9.2     2023-02-10 [?] CRAN (R 4.2.0)
#>  P magrittr               2.0.3     2022-03-30 [?] CRAN (R 4.2.0)
#>  P Matrix                 1.5-4.1   2023-05-18 [?] CRAN (R 4.2.0)
#>  P MatrixGenerics       * 1.10.0    2022-11-01 [?] Bioconductor
#>  P matrixStats          * 0.63.0    2022-11-18 [?] CRAN (R 4.2.0)
#>  P metapod                1.6.0     2022-11-01 [?] Bioconductor
#>  P mgcv                   1.8-42    2023-03-02 [?] CRAN (R 4.2.0)
#>  P mime                   0.12      2021-09-28 [?] CRAN (R 4.2.0)
#>  P munsell                0.5.0     2018-06-12 [?] CRAN (R 4.2.0)
#>  P nlme                   3.1-162   2023-01-31 [?] CRAN (R 4.2.0)
#>  P pillar                 1.9.0     2023-03-22 [?] CRAN (R 4.2.0)
#>  P pkgconfig              2.0.3     2019-09-22 [?] CRAN (R 4.2.0)
#>  P plyr                   1.8.8     2022-11-11 [?] CRAN (R 4.2.0)
#>  P promises               1.2.0.1   2021-02-11 [?] CRAN (R 4.2.0)
#>  P purrr                * 1.0.1     2023-01-10 [?] CRAN (R 4.2.0)
#>  P R6                     2.5.1     2021-08-19 [?] CRAN (R 4.2.0)
#>  P Rcpp                   1.0.10    2023-01-22 [?] CRAN (R 4.2.0)
#>  P RCurl                  1.98-1.12 2023-03-27 [?] CRAN (R 4.2.0)
#>  P readr                * 2.1.4     2023-02-10 [?] CRAN (R 4.2.0)
#>    renv                   0.17.3    2023-04-06 [1] CRAN (R 4.2.0)
#>  P reshape2               1.4.4     2020-04-09 [?] CRAN (R 4.2.0)
#>  P rlang                  1.1.1     2023-04-28 [?] CRAN (R 4.2.0)
#>  P rmarkdown              2.21      2023-03-26 [?] CRAN (R 4.2.0)
#>  P ROCR                   1.0-11    2020-05-02 [?] CRAN (R 4.2.0)
#>  P rprojroot              2.0.3     2022-04-02 [?] CRAN (R 4.2.0)
#>  P rstudioapi             0.14      2022-08-22 [?] CRAN (R 4.2.0)
#>  P rsvd                   1.0.5     2021-04-16 [?] CRAN (R 4.2.0)
#>  P rvest                  1.0.3     2022-08-19 [?] CRAN (R 4.2.0)
#>  P S4Vectors            * 0.36.2    2023-02-26 [?] Bioconductor
#>  P sass                   0.4.6     2023-05-03 [?] CRAN (R 4.2.0)
#>  P ScaledMatrix           1.6.0     2022-11-01 [?] Bioconductor
#>  P scales                 1.2.1     2022-08-20 [?] CRAN (R 4.2.0)
#>  P scater                 1.26.1    2022-11-13 [?] Bioconductor
#>  P scran                  1.26.2    2023-01-19 [?] Bioconductor
#>  P scuttle              * 1.8.4     2023-01-19 [?] Bioconductor
#>  P sessioninfo            1.2.2     2021-12-06 [?] CRAN (R 4.2.0)
#>  P shiny                  1.7.4     2022-12-15 [?] CRAN (R 4.2.0)
#>  P shinyBS                0.61.1    2022-04-17 [?] CRAN (R 4.2.0)
#>  P shinydashboard         0.7.2     2021-09-30 [?] CRAN (R 4.2.0)
#>  P SingleCellExperiment * 1.20.1    2023-03-17 [?] Bioconductor
#>  P sparseMatrixStats      1.10.0    2022-11-01 [?] Bioconductor
#>  P statmod                1.5.0     2023-01-06 [?] CRAN (R 4.2.0)
#>  P stringi                1.7.12    2023-01-11 [?] CRAN (R 4.2.0)
#>  P stringr              * 1.5.0     2022-12-02 [?] CRAN (R 4.2.0)
#>  P SummarizedExperiment * 1.28.0    2022-11-01 [?] Bioconductor
#>  P svglite                2.1.1     2023-01-10 [?] CRAN (R 4.2.0)
#>  P systemfonts            1.0.4     2022-02-11 [?] CRAN (R 4.2.0)
#>  P tibble               * 3.2.1     2023-03-20 [?] CRAN (R 4.2.0)
#>  P tidyr                * 1.3.0     2023-01-24 [?] CRAN (R 4.2.0)
#>  P tidyselect             1.2.0     2022-10-10 [?] CRAN (R 4.2.0)
#>  P tidyverse            * 2.0.0     2023-02-22 [?] CRAN (R 4.2.0)
#>  P timechange             0.2.0     2023-01-11 [?] CRAN (R 4.2.0)
#>  P tzdb                   0.4.0     2023-05-12 [?] CRAN (R 4.2.0)
#>  P UpSetR                 1.4.0     2019-05-22 [?] CRAN (R 4.2.0)
#>  P utf8                   1.2.3     2023-01-31 [?] CRAN (R 4.2.0)
#>  P vctrs                  0.6.2     2023-04-19 [?] CRAN (R 4.2.0)
#>  P vipor                  0.4.5     2017-03-22 [?] CRAN (R 4.2.0)
#>  P viridis                0.6.3     2023-05-03 [?] CRAN (R 4.2.0)
#>  P viridisLite            0.4.2     2023-05-02 [?] CRAN (R 4.2.0)
#>  P webshot                0.5.4     2022-09-26 [?] CRAN (R 4.2.0)
#>  P withr                  2.5.0     2022-03-03 [?] CRAN (R 4.2.0)
#>  P xfun                   0.39      2023-04-20 [?] CRAN (R 4.2.0)
#>  P xml2                   1.3.4     2023-04-27 [?] CRAN (R 4.2.0)
#>  P xtable                 1.8-4     2019-04-21 [?] CRAN (R 4.2.0)
#>  P XVector                0.38.0    2022-11-01 [?] Bioconductor
#>  P yaml                   2.3.7     2023-01-23 [?] CRAN (R 4.2.0)
#>  P zlibbioc               1.44.0    2022-11-01 [?] Bioconductor
#> 
#>  [1] /kyukon/data/gent/vo/000/gvo00063/Jeroen/DD_project_v2/benchmarks/lupus-n_patients/renv/library/R-4.2/x86_64-pc-linux-gnu
#>  [2] /kyukon/home/gent/460/vsc46052/.cache/R/renv/sandbox/R-4.2/x86_64-pc-linux-gnu/4df86545
#> 
#>  P ── Loaded and on-disk path mismatch.
#> 
#> ──────────────────────────────────────────────────────────────────────────────

Lupus: influence of number of patients - mock results

Results for celltypes = B_mem

Milan Malfait & Jeroen Gilis

June 02, 2023

Setup

Load results

Load SCE objects

Data overview

n_patients = 10

replicate_1

replicate_2

replicate_3

replicate_4

replicate_5

n_patients = 20

replicate_1

replicate_2

replicate_3

replicate_4

replicate_5

n_patients = 30

replicate_1

replicate_2

replicate_3

replicate_4

replicate_5

Extract results of interest

Runtimes

P-values

Visualize results

Run times

Number of significant genes for different FDR levels

n_patients_10

n_patients_20

n_patients_30

P-value distributions

n_patients_10

n_patients_20

n_patients_30

P-value distributions stratified by expression level

n_patients_10

n_patients_20

n_patients_30

Session info

Results for `celltypes =` B_mem

`n_patients = 10`

`n_patients = 20`

`n_patients = 30`